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The reduced dynamics formalism has recently emerged as a powerful tool to study the dynamics 
of nonequilibrium quantum impurity models in strongly correlated regimes. Examples include the 
nonequilibrium Anderson impurity model near the Kondo crossover temperature and the nonequi- 
librium Holstein model, for which the formalism provides an accurate description of the reduced 
density matrix of the system for a wide range of timescales. In this work, we generalize the formalism 
to allow for non-system, observables such as the current between the impurity and leads. We show 
that the equation of motion for the reduced observable of interest can be closed with the equation 
of motion for the reduced density matrix and demonstrate the new formalism for a generic resonant 
level model. 



I. INTRODUCTION 

The study of open quantum impurity models, where 
the coupling of a small system to multiple baths drives it 
permanently away from the possibility of an equilibrium 
state, is an active and rapidly progressing field of re- 
search. It has recently become possible to make quantita- 
tive statements about experimentally measurable trans- 
port properties in certain cases ;i~' however, in general 
many unresolved issues remain. For example, the na- 
ture of the charge and spin dynamics within the Kondo 
regime of a quantum dot driven out of equilibrium is cur- 
rently under investigation,^ and basic questions regarding 
hysteresis and bi-stability in systems governed by strong 
electron-phonon couplings remain under debate. ISU^ A 
major theoretical challenge lies in the need to provide an 
accurate account of time propagation of such open quan- 
tum systems, starting from some known initial state and 
proceeding all the way to an unknown steady-state. 

Numerically exact methods play a particularly impor- 
tant role in the quest to obtain a reliable, unbiased de- 
scription of nonequilibrium phenomena. Several differ- 
ent types of brute-force approaches developed in recent 
years have been applied to open nonequilibrium quan- 
tum systems. These include the time-dependent numer- 
ical renormalization-group technique ,112113 jterative^^HiSl 
or stochastie^^^ dia gramm atic methods, and wavefunc- 
tion based approaches .1221231 ■\Yhile the application of these 
approaches to the the nonequilibrium Holstein, the An- 
derson impurity, and the spin-fermion models has been 
very fruitful, they are still restricted to a relatively small 
range of parameters, typically characterized by a rapid 
decay to steady-state. Situations or observables exhibit- 
ing slow dynamics are inaccessible by these brute-force 
methods. 



An alternative approach recently proposed by Cohen 
and RabanPS is based on a combination of a brute- 
force impurity solver (one of the above) with a general- 
ized quantum master equation (GQME). The Nakajima- 
Zwanzig-MorP^^^ formalism was used to derive an ex- 
act equation of motion for the reduced density matrix 
of the system, which includes a memory kernel giving 
rise to non-Markovian effects. This kernel, along with 
some information regarding the initial conditions, deter- 
mines the dynamics of the system and contains all infor- 
mation about the time dependence of single-time system 
observables and their steady-state values. In many sit- 
uations of interest, in particular when the bandwidth of 
the baths is large compared to other energy scales in the 
proble m, the m emory kernel is expected to decay rapidly 
to zero.'^'23E3 Thus, one can safely truncate the memory 
kernel at a finite time, performing a "cutoff approxima- 
tion". Brute-force impurity solvers limited to short times 
are well suited for the kernel's numerical evaluation up 
to the cutoff time, and once the memory kernel has been 
obtained, the GQME is exact and tractable at all times. 

The GQME formalism has rece ntly been combined 
with the Bold impurity solvei'2i'29l ^q uncover the spin 
dynamics near the Kondo crossover temperatur&B and 
with the multilayer multi-configuration time-dependent 
Hartree method to reveal the nature of bi-stability in 
systems with electron-phonon couplings .1221 Despite the 
open nature of the systems studied in these works, trans- 
port properties were not addressed; this is due to one 
of the formalism's main limitations, in that observables 
outside the impurity part of the Hilbert space are not 
accessible, and only system observables such as the dot's 
population or magnetization can be calculated. On the 
other hand, perturbative expressions for transport prop- 
erties in terms of vertex functions have been derived and 



evaluated b efore in approximate methodologies built on 
the GQMEJSSEB and it seems reasonable to expect that 
a general exact formulation in the spirit of Ref|231 should 
exist. 

In this paper we extend the GQME formalism (re- 
viewed in Section In]) to describe non-system observables. 
This allows for comparison of predictions made by the 
GQME with a much wider variety of experimental ob- 
servables, of which an important example (worked out 
in detail here) is the current. Equally important, it fa- 
cilitates access to the spectral functions by way of mea- 
suring the current to an infinitesimally coupl ed a uxiliary 
bath|22H2ll xhe key idea discussed in Section III is based 



on deriving a reduced equation of motion for the observ- 
able of interest, which can then be expressed in terms of 
the reduced density matrix. This leads to an introduc- 
tion of an additional, observable-specific memory kernel 
with properties qualitatively similar to those of the mem- 
ory kernel appearing in the standard GQME. Section IV 
is devoted to expressing the steady-state properties in 
terms of the memory kernels alone, while in Section [V] 
we show how the projected quantities appearing in the 
GQME can be translated into the language of ordinary 
observables expressed as second-quantized operators, us- 
ing a noninteracting model as an illustrative example. In 
Section |VT] we present several test cases and examples for 
the non-interacting case, where the properties of memory 
kernels can be explored without the need for technically 
complicated numerical solvers. Finally, a summary is 
given in Section |VlT] 



II. PROJECTED DYNAMICS FOR SYSTEM 
OBSERVABLES 

We will begin by reviewing the derivation of exact pro- 
jected (or reduced) equations of motion,^- and the pro- 
cess of going from projected to unprojected dynamics. ES 
These details are provided here in a self-contained man- 
ner because they will be of particular importance later, 
when we discuss how the process can be generalized. 
Consider an operator Hilbert space H = S(E)B composed 
of two subspaces S and B, which we will call the system 
and bath subspaces. We are interested in a Hamiltonian 
of the form 



H^Hs + Hb+V, 



(1) 



where Hs £ S is the system or impurity Hamiltonian, 
Hb E B is the bath Hamiltonian and V E H, V ^ S, B, 
is the coupling Hamiltonian. Generally, the motivation 
for employing such a description is to describe a small, 
strongly interacting impurity coupled to large noninter- 
acting baths, but we need make no further assumptions 
at this stage. We can now define a projection operator P 
onto the S subspace by tracing out the bath degrees of 
freedom, in the process also defining its complementary 



operator Q: 



Here ps = e 



P = pb^^b, 
Q^l-P. 



(2) 
(3) 



-PHe 



/Trse" 



-I3He 



We also define ps E S to 



be the initial impurity density matrix, and po = Pb (^ PS 
the initial full density matrix. The expectation value of 
a system operator ^ S 5 is given by 



{A (t))^ Tip it) A 

= Trs[(TrBp(t))A] 
= Trs{a(i)A}, 



(4) 
(5) 
(6) 



and the reduced density matrix a (t) = TtbP (t) contains 
information about all single-time properties of system ob- 
servables. This object has a lower dimensionality than 
that of p, and it would thus be economical to describe its 
equations of motion without referring to the system as 
a whole. This is the basic idea behind reduced quantum 
dynamics. 

To proceed, one considers the Liouville-von Neumann 
equation, which governs the dynamics of the full density 
matrix: 



'^dtP 



[H,p] = Cp. 



(7) 



The Liouvillian superoperator C denotes performing a 
commutation with the Hamiltonian, such that CA = 
[H,A]. We also define CsA = [Hs,A], CyA = [Hv,A] 
and CbA = [Hb,A]. Applying each of the projection 
operators from the left and using 1 — P + Q within the 
commutator gives: 



ih--Pp = P[H,{P + Q)p] 
at 

ih^Qp = Q[H,iP + Q)p] 
at 



(8) 
(9) 



Eq. (l9| has the formal solution 

Qp = e-^^^'Qpo 

f dre-'Q^^QCpBait-T), 
■ Jo 






(10) 



which can be inserted into Eq. (Is]) to obtain the 
Nakajima-Zwanzig-Mori equation (NZME)BM2Z1 



iha [t) = Cscf it) +i9(t)- - / dTK (r) a {t - t) /ll) 

^ Jo 

(12) 
(13) 



K(t) = TYB{Cve-''^'^^QCpB}, 
^(t) = TrB{£ye-HO^*Qpo}. 



Let us take a moment to examine the important re- 
lation of Eq. (111. It has the form of an operator lin- 



ear Volterra integro-differential equation of the second 
kind. As we have made no approximations, it is exact; 
yet it contains only operators and superoperators within 
the low-dimensional system space. The time derivative 



of the reduced density matrix a is given by the sum of 
three contributions: the first term (£scr (i)) describes the 
exact evolution of the system if the couphng to the bath 
were set to zero. The second term (^ (i)) expresses ini- 
tial correlations between the system and bath, and it is 
easy to verify from its definition in Eq. ( |13[ ) that it equals 
zero for the factorized initial conditions po = Pb ® Ps- 
The last term includes the memory kernel {k{t)), and 
depends on the complete history of a (t) at earlier times. 
The appearance of this non-Markovian term is the price 
of going to reduced dynamics, and to make headway with 
the NZME one must begin by evaluating k (t) . 

The definition of k (i) in Eq. ( 12 1 includes the trou- 
blesome component e^^Q-^'^. To understand why it is 
troubling, consider the following: it is easy to show that 
the superoperator e~i^'^ evolves the density matrix with 
respect to the Hamiltonian, thus simply expressing our 
familiar notion of dynamics: 



Eqs. (18 1 and (21 1 are superoperator linear Volterra 



-iiCr , 



AHt 



pe 



„-iH7 



(14) 



The modified operator e^f^'^^'^ , however, contains a pro- 
jection operator within the exponent, and so does some- 
thing else entirely — something which turns out to be sub- 
stantially harder to understand or calculate. Our next 
step is therefore to get rid of these inconvenient projected 
dynamics. While several ways to go about this task ex- 
ist, we will limit the discussion to a particular method 
suggested by Zhang et al)^ 

Consider the function ^ (t) of Eq. (13 I. By applying 
the identity 



g-iQCt ^f,-ict 



- f dTe-i^'^'-'^^PCe-iQ^-' (15) 
T-Jo 



to its definition, we can obtain: 



it) = Trs { 



Cye " 



^^*Qpo} 



(16) 



TYB<Cve-^^'Qpo 



j_ I dr£ye-*^(*-")p£e-*^^"Qpo|(17) 



where 



= S(0-$(t)a(0) 

+ "- [ dT$ (t - r) i9 (t) , 

" JO 

E{t)=TTB{Cve-i^'po], 
<i>{t)^TTB[Cve-i^'pB}. 



(18) 

(19) 
(20) 



Applying the same identity Eq. (15 1 to Eq. (12 I yields 



K (t) = ih<^ (t) - $ (r) £5 + T / dr$ {t - t) k (t) . (21) 



integral equations of the second kind, with both the in- 
homogeneous contributions and the kernels determined 
by the combination of Eq. ( 19 I and ( 20 1 and the form of 



the system Liouvillian operator. Like the NZME, Eq. [TT] 
they consist of objects which inhabit the low-dimensional 
impurity subspace — however, they have a higher dimen- 
sionality due to their superoperator nature (if a can be 
represented by an A'^ x iV matrix, then •& and k are 
iV^ X A^). Their importance lies in the fact that $ and 
S, which are propagated by the full Hamiltonian with 
normal dynamics, can be written in terms of physical 
observables; this means they can be evaluated with a va- 
riety of computational methods, and then used to solve 
Eqs. (18 I and 21 numerically. 



We now have the necessary description to introduce 
the cutoff approximation: if we have some way of evalu- 
ating K {t) up to some finite time, it is sometimes possible 
to make an ansatz about later times. Importantly, if the 
memory has decayed to zero to within a numerical accu- 
racy over a finite time, one can proof that it will remain 
zero at all later times. One then solves the NZME with 
this cutoff memory kernel to obtain an approximate value 
for (7 (t): however, if a{t) can be converged in the cutoff 
time to within the desired accuracy, the entire procedure 
is numerically exact. 



III. GENERALIZED PROJECTED DYNAMICS 

FOR NON-SYSTEM OBSERVABLES 

The time-dependent electronic current flowing through 
an impurity does not have a single definition, as it de- 
pends in general on the topology of the surface through 
which electronic flow is measured. In steady state popu- 
lations must be constant, and the current must therefore 
become independent of this definition (if it is unique); 
however, the definition itself remains arbitrary. If the 
impurity model in question is, for instance, given by 
a chain Hamiltonian, current may be measured at any 
point along the chain. However, in models where not all 
current must fiow between impurity sites, it is customary 
to measure currents at the junction between the impu- 
rity and one of the leads (baths). In a setup involving two 
Fermionic baths held at different chemical potentials, of- 
ten referred to as the left (L) and right (i?) leads, one 
is therefore interested in quantities such as the so-called 
"left current" (or alternatively the "right current"): 



It 



dt 



-eN, 



qeL 



[H, alag] 



(22) 



Here the aj, and a/, are creation and destruction operators 
in the left lead subspace L C B. The current operator 
will therefore in general not be a member of the impurity 
subspace S, and cannot be obtained from Eq. ([6| along 
with knowledge of the reduced density matrix a(t). It 
should be mentioned briefly that one simple way of deal- 



ing with this issue is to define an effective Hamiltonian 
and a repartitioning of "H in such a way that the cur- 
rent can be measured within the system; however, this 
will invariably raise the dimensionality of S, which may 
complicate the problem beyond the applicability of many 
numerical methods. 

In order to allow for the calculation of non-system ob- 
servables, we proceed by deriving a Nakajima-Zwanzig- 
Mori-like equation for the expectation value of a general 
operator I. While / implies that we are interested in the 
current, nothing in this section is limited to that specific 
case. The ideas to follow could work equally well for any 
operator, but for the current one might expect them to 
converge quickly with a cutoff time, since it is expected 
to be determined largely by quantities local to the dot. 

We start from the equation of motion for / in the 
Schrodinger picture, where p has a time dependence but 
/ does not: 



th-ip = i[H,p] 



(23) 



Applying the projection operators using the definitions 
and procedure of the previous section yields: 



ih—PIp = PICPp + PICQp, (24) 

ih—QIp = QICPp + QICQp. (25) 

In addition to these, we still have equations (l8| and ^ 
for the density matrix, which may be written in the form: 



ih—Pp = PCPp + PLQp, 
at 

ih^Qp = QCPp + QCQp. 
at 

The latter has the formal solution: 



(26) 

(27) 



Qp^e~iQ^'Qpo 

/ dre-'Q^(*-")Q£pBa(T). (28) 
Jo 



h 



Putting this expression together with Eq. ( 24 ) and defin- 
ing L (t) = Trslp, we obtain 



ihl (t) ^ TtbICpbct (t) + TiBlCe-i^'^^Qpa 



(29) 



drTrs {lCe-''^^'^'-^^QCpB\ o (i 



or 



ihi (t) = C,a (t) +^,{t)-- / dr K, {t^T)a (r) . (30) 



The terms of this equation appear similar to those of 
the NZME in Eq. (Ill, though it is a solution in closed 



L (t) appears only on the left hand side. In writing it we 
have defined: 

C, = Ttb{ICpb}, (31) 

d,{t)=TTB{lCe-iQ^'Qpo}, (32) 

K, (t) = Trs {iCe-'Q'^'QCpB } • (33) 

The initial correlation term "fJ^ is once again zero for 
uncorrelated initial conditions and this time we will re- 
move it for the sake of brevity. As in the formalism for cr, 
in order to phrase everything in terms of quantities with 
unprojected dynamics we now once again perform the 
Zhang-Ka-Geva transformation^ on the current mem- 
ory kernel k^. Applying the identity (151 allows us to 
write 



K,{t)^TTB{lCe~i'^'QCpB} 



drTrR|//:e-s^('-^)p£e 



= ^;^$,, (t) - $, (t) Cs 

rt 

dr$t {t~ t) k{t) 



(34) 
(35) 



or 



K, (t) - ^;^$, (t) - $,, (i) Cs 

dr$,, {t- t)k (t) 



(36) 



Once again, this is a closed form solution rather than an 
integral equation. Its inputs are the same k (r) defined 
in Eq. (12 I, as well as the new quantity 



<^,{t)=Ti-B{l£e-i'^'pB} 



(37) 



Eqs. (36 1 and (30 1 amount to a generalization of the 



Nakajima-Zwanzig-Mori formalism to non-system oper- 
ators. The structure of these equations is reminiscent of 
the structure of the corresponding equations in the orig- 
inal theory, on which the extension relies — and yet they 
are simpler in a certain sense, as they are closed form so- 
lutions up to quadrature rather than integro-differential 
or integral equations. In addition to a (t) and k (t), which 
can be obtained from the original theory, the extended 
formalism relies on a new input, $,, (t), which is defined 
in terms of regular (rather than projected) time prop- 
agation and must be calculated explicitly. Once $t (t) 
is available one can solve Eq. ( |36[ ) to obtain k, (i), and 
then solve Eq. (30 1 to obtain t (i), a system-space oper- 
ator which can be traced over to obtain the expectation 
value of the operator /. 

IV. STEADY STATE 

If we wish to examine the i — > oo limit of a (i), it is 
more convenient to define the Laplace transform 



form rather than an integro-differential equation, since 



a(z) 



e-'-*a (t) dt 



(38) 



When applied to Eq. (11), this yields: 



ih [za {z) - a (0)] = Lgd {z) + d {z) - -k {z) a (^)(39) 



a{z)^ 



a{Q) + i-J{z) 



- -C 



rk (z) 



(40) 



Using the final value theorem a {po) — lini2_>.o za (2), we 
can obtain an expression for a at long times: 



cr(t 



iha (0) + 1? (z) 



lim 



(41) 



We can also obtain a stationary-state equation by consid- 
ering a time independent solution cr (t — >■ 00) to Eq. (Ill 



such that we can set the time derivative to zero and take 
a outside the integral before taking the Laplace trans- 
form. If we also assume that the initial correlations are 
either zero to begin with or die out at infinite time, this 
gives: 



Ls--^k{z 



iQ) ]a{t^oo)=Q. 



(42) 



This last equation is of particular interest, because it 
allows us to go from the memory kernel and system Li- 
ouvillian directly to the steady state properties of the 
reduced density matrix, without passing through the dy- 
namics and without any reference to the initial state or 
correlations of the system. This is very useful when we 
are interested in general questions regarding the steady 
state, such as that of its existence or uniqueness. When 
applying the cutoff approximation, k (z — >■ iO) must be 
calculated to sufficient accuracy that the steady state 
density matrix converges. 

It is natural to attempt deriving a similar expression 
for the current directly at steady state. One way of going 



about this task is to begin with Eq. ( 30 1 and take the 
Laplace transform: 



ih zl (z) - i (0) = £,<T (z) + d, (z) 



--«,(z)a(z). (43) 

Extracting zl (z) and using the final value theorem then 
gives 

i{1.^^)^\mi\\L,-\k,{z)\a{z) (44) 

This suggests that in order for a steady state to exist, we 
must have 



lim £t — —ki (z) ~ z. 



The constant of proportionality (itself a superoperator) 
determines the value of the current at steady state. In 
the case of impurity observables, it is sufficient to know 
the zero frequency component of the memory kernel in 
order to obtain the steady-state value. Here, however, 
one must also know something about the low-frequency 
properties (or linear frequency response) of the current 
memory kernel, k^ (z). 



V. EXPRESSING THE KERNELS IN 
SECOND-QUANTIZED FORM 

Everything up to this point has been independent of 
the details of any particular model, requiring only that 
a partitioning between the impurity and bath part be 
made. In order to illustrate the process of using the for- 
malism presented above in a particular model, we will 
continue by way of the simplest possible example: that 
of a noninteracting junction (the formalism is not limited 
to this case3l2iEI!) . This model, often called the resonant 
level model, is defined by the Hamiltonian 



H = Hs + Hl + V, 


(47) 


Hs = edU, 


(48) 


Hl = ^Egflja,, 


(49) 



y = ^tgdaj+t;a,d^ 



(50) 



A complete definition must include the e^and tq, and in 
this case all necessary information is contained in the lead 
coupling function 



2n}^\mUiu^- 



(51) 



The first step in the calculation is the evaluation 
of the system Liovillian. It is convenient to work in 
the Hubbard representation for operators in the impu- 
rity subspace: an operator A G S can be written as 
A = '^^^ ciij \i) (j'l, where the indices i and j can take 
on the values of states in the impurity subspace — in this 
case and 1 for unoccupied and occupied, respectively. 
This superoperator simply performs a commutation with 
the system Hamiltonian, and using Eq. (48) to insert the 



explicit form of Hs yields an expressions for Cs in matrix 
(or tetradic) form: 



[-^«]^J^fc^ = Tr 



s{i\z){j\)^Cs\k){l\} (52) 

1 
Y,{m\{\^){J\)Us\l){l\,\k}{l\]\m) (53) 

= £ [SjioSiki + Sijkii - SjiiSik] ■ (54) 

(4dj When no bath is present the model is reduced to a two- 



level system, and Eq. (|11| gives the expected result: 

- o-fc/ (55) 



^^ ki 



= £ (^jif^joo-io - 5iQ6jiaQi.) (56) 

That is. ofF-diagonal density matrix elements oscillate 
with a frequency | while diagonal elements remain sta- 
tionary. 

Next, we need to evaluate the memory kernel. This 
requires the evaluation of the superoperator 

$(t)A = TrB{£ve-^^*PB^} (57) 

= Trs ive-i"^pBAeT^"^ 

-e-^^^psAe^^^v] , (58) 

which can also be represented in matrix form: 

^..,Mit)^TTs{i\^){J\)^c^it)\k){l\} (59) 

= A-A\ (60) 

with 

A ^ Ti's {\j) {i\ Ttb {Ve-i'^'pB \k) {l\ e*^*}} ,(61) 

A' EE Trs {\j) {^\ Tis {e-i"'pB \k) {l\ ei"'v}} .(62) 

Consider the term A. Let us perform the trace over the 
impurity space and take the Hubbard operators to second 
quantized form: 

A = TvB {pB {l\ e^""* 

[SiiSjid'^d + d^oSjodd'^ + S^iSjod + SioSjid^) 
Ve-i"'\k)} (63) 

= TtbIpb {l\ 

{SaSjid^t)d{t) + S,oS,od{t)d^t) 

+SaSjod{t)+SMdjid^t)) 

Vit)\k)]. (64) 

In the final step, the operators were given their full time 
dependence in the Heisenberg picture. Using V (i) = 
J2n ^qd (t) flg it)+t*qaq (t) d^ (t) and the fact that all pairs 
of dot and lead operators maintain normal commutation 
relations when taken at the same times, one can now 
show that 

A = J2 ^^Trs {pB {l\ S^Sjod (t) aj (t) \k)} (65) 

+ ^ t.Trs {pB {l\ S.oS.id^ {t)d{t) 4 it) \k)} 
q 

- E K^^B {pB {l\ 6,i5,^d^ [t) a, {t) \k)] 
q 

-Y^t^TiB {pB {l\S,i6jod{t)d^ {t)ag{t) \k)} . 



Similarly, 

A' ^Y^tqTvB {pB {l\ Sa6jid{t) al (t) \k)} (66) 

q 

- ^t^Trs {pB {l\ S,o6jid{t) S (t) al {t) \k)} 

q 

- E^^Trs {pB {l\ S^oSjod^ {t) a, (t) |fc)} 

q 

+ ^t;TrB {pB {l\S,iSjod^ {t)d{t)ag (t) \k)} . 



Putting the expressions for A and A' into their defining 
equation then yields 






(67) 



with 



VM = Trs I ^ tqPB {l\ d {t) a\ (t) \k) I , (68) 
i^ki=TTB\Y.tqPB{l\al{t)\k)\. 



(69) 



The if elements have a rather simple physical interpreta- 
tion: they are directly proportional to the time derivative 
of the total population on the dot. The ijj elements are 
harder to interpret in such a manner. 

The equations therefore collapse to a simple form, 
phrased in terms of the system-space matrix elements 
ipki S'lid ipki of normal second quantization operators 
propagated under the influence of the full Hamiltonian. 
Some further simplification can be made by considering 



Eq. ( 68 1 as a function of time when we go to the interac- 
tion picture under Hq = Hg + Hb '■ 



(70) 



^H=TrBlY.^qPB{l\e'"'dale~^'''\k)\ 

= Y, igTrs 

q 

{pB{l\UHt)dHAt)a^Ho.q(t)U{t)\k)}. (71) 

It is easy to verify that the time dependence of d and a' 
in the interaction-picture is described by a simple oscil- 
lation, and that the U and W are sums over products 
of terms containing either a^d or d^aq in the interaction 
picture. The trace over the bath may then be performed 
at time zero, throwing out all terms which do not have 
the same number of Uq and aj operators. Yet, from the 
argument we have just stated, such terms will also have 
the same number of d and d^ operators, and (pki must be 
zero unless k = I. For similar considerations, ipj^i is zero 



unless k ^ I and we have 



First, consider the current LiouvilHan operator: 



Vki = Trs J 51 ^PPB (P\ d it) 4 (t) \k) \ Ski, (72) 
i^ki = Trs I 5] t,pB {l\ al (t) \k) I (1 - Ski) ■ (73) 



Considering the role of ip and ^jj in Eq. (671, one can 



see that $ contains terms which couple the populations 
and terms which couple the coherences; it contains no 
terms which couple the populations to the coherences. 
Examining Eq. (211, one realizes that terms with no in- 



homogeneous contribution must identically vanish, such 
that 



^ij^kl 






(74) 



otherwise. 



Similarly, Eq. (Ill, along with the Liouvillian (54 1, im 



mediately leads us to the conclusion that the diagonal 
elements of a form one coupled block within the formal- 
ism, while the off-diagonal elements of a form a second 
block: in other words, within the resonant level model, 
the reduced dynamics of the diagonal elements (the pop- 
ulations) are decoupled from those of the off-diagonal el- 
ements (the coherences). 

Among other things, this implies that if we are inter- 
ested only in the populations we do not need to calcu- 
late the ipki, and vice- versa for the coherences. Since 
the populations are also unaffected by the Liouvillian, 
and assuming factorized initial conditions, the equation 
of motion turns from a superoperator equation into a 
matrix equation for the population vector an: 



ih — CTr, 

at 



it) = - 



dr y K; 

j 



"jj 



(T)a„(t-r). (75) 



Interestingly, this analytic block structure conclusion 
holds in the generalized Holstein modeP^ as well (though 
this will not be shown here), and continues to hold for 
the Anderson modeP* even in the presence of a magnetic 
field.A 

We continue to examine the memory formalism for the 
non-system current operator in the noninteracting case. 
We will define a simplified left current operator as 



J^tqdal, 

qeL 



(76) 



such that I (t) = TiB \ Ip (t) r and the physical current 

is 

(/ (t)) = -2^3 (/ (t)) = -2^5Tr5 {im . (77) 



Trs liCpBJ = Tr_B llCpEJ 



(78) 



Tib \ J2 hda\CvPB ) ■ (79) 

The term containing Cb can be dropped because 
CbPb-A. = {CbPb)A = for any system variable A, 
while the term containing Cg is zero because it must 
contain a trace over an odd number of lead creation or 
destruction operators. It is then simple to show by the 
same procedure from which the system Liouvillian was 
derived that 



[A],,- ,„, = Trs {(|z) (j|)^ Trs [iCpb] \m) {l\] (80) 

= ^|t,|'Tr{(|*) {j\)Ud^a\a,pB \m) (?|}(81) 

q£L 

~Y.\t,\'i:T[{\z){j\)Ua\pB\m){l\a,d^] 

qeL 

= Y.U^^TB{{i\ddU,PB\Tn) (Z|j)} (82) 



- X! I*«l^ ^^^ i^*! do^tP^ 1'^^ ^'1 "9*^^ l-^'H ' 

q£L 

where f„ = stt — tt^ is the Fermi-Dirac distribution. 

In the above equation, we have used the factorized initial 
conditions by assuming an equilibrium Fermi distribution 
at i = in the baths (it is worth noting that this distri- 
bution is allowed to evolve freely under the influence of 
the full Hamiltonian in the reduced dynamics formalism, 
yet the full details of bath dynamics are no longer acces- 
sible from the information stored in a (<)). With this, it 
is straightforward to show that 



[-C.],,,™, = S,o5moSi,l^< (0) - (5,o5,o<5mi(5aA> (0) (83) 
where the hybridization functions 

A< (t) ^ I du e-* Y. 1^*1' f<i^ (^ - ^?) ' (84) 

•' qeL 

A> (t) - / dc. e'"* J2 l^?l' (1 - /«) ^ ('^ - ^^9) ' (85) 
"^ qeL 



can easily be evaluated in terms of the frequency-space 
coupling density given in Eq. (51). 



The final object we need to evaluate is <l>t. Since no 
new conceptual issues arise here as compared with the 
calculation performed for $, we will simply write down 



the final answer. With the definitions 



^f.l 



E 

qeL.q' 



(86) 



-Y.[\U\' {dHt)d{t)) +t,e,{d{t)al{t))) , 

q£L 

V.,2^Y.^q{d{t)a\{t)), (87) 

qeL 

^-,,,1= E kt*^'{dHt)al{t)a,,{t)), (88) 

q£L,q' 

^.,2= J2 Uh'{da\,{t)a\{t)), (89) 

q£L,q' 

<i>t (i) takes the simple form: 

+<5.o'5,i (V'M - ^-.,2)} |A:) . (90) 

The inherent asymmetry of the expression above is due 
to the asymmetric definition of /. As before, it is easy 
to show that the block structure is such that the current 
can be determined without reference to the off-diagonal 
elements of either (j{t) or t(i). Furthermore, it is straight- 
forward to show that for the resonant level model (and 

for the Holstein model) [<i>t]gQ qq (t) — ff(0) for an ini- 
tially empty dot and [$(,]oo n (i) = ft (1) f^^' ^^^ initially 
occupied dot. These are useful relations as they provide 
an alternative way of computing k^ (t) directly from the 
left current (at short times) , without the need to evaluate 

VI. RESULTS 

The physics of the resonant level model are generally 
well known, yet the literature has seen little exploration 
of the properties of the memory kernel in this model, and 
of course none of the current memory kernel which has 
been introduced here. We therefore present some results 
below which we expect to be of interest to the field, as 
they provide insight into those aspects of the problem 
which do not rely on interaction. In order to restrict the 
parameter space explored, we will discuss the symmetric 
case in which e = with a bias voltage applied symmet- 
rically such that V = 2/i/^ = —2/1^ (from here on we 
set h = e = I). The lead coupling densities are taken 
to be Tl.r{^) = 7 7 — „ n ] 7 ?rTV- We will 

:i\ / (l+e"('^-'V)j(l+e"(-"-"c)j 

limit our attention only to the diagonal elements of the 
reduced density matrix and the corresponding element of 
the memory kernel; these elements are completely decou- 
pled from the off-diagonal coherences, as discussed above, 
and therefore no approximation ensues from this. 

All the results presented in this section are exact and 
have been calculated by a direct solution of the full equa- 
tion of motion of the complete density matrix, a tech- 
nique which relies on the quadratic form of the Hamilto- 
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Figure 1. An element of the memory kernel k of the fully sym- 
metric, resonant level model at a range of band parameters. 
Due to the symmetry, the results shown here are independent 
of both temperature and voltage. Panel A through C show the 
effect of softening the band edge by varying u, while within 
each separate panel the effect of varying the bandwidth (which 
is approximately twice the cutoff frequency Qc) is illustrated. 
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Figure 2. The two nonzero elements of k^ for the left current 
are shown in panels A and B, with Tf = 10, r/3 = 1 and 
V = 4r. In each panel, the time dependence of the k,, element 
is shown at a range of band widths. 



nian and is therefore applicable only to the noninteract- 
ing case. In general, making similar progress for inter- 
acting s ystems req uires a numerical solver of one type or 
anotherPHMIMzr 

We begin with a discussion of the behavior of the mem- 
ory kernel. The symmetrical parameters we have chosen 
are of particular interest because in the absence of inter- 
action both a (t) and k (t) are completely independent 
of both the temperature and voltage. In addition, all 
nonzero matrix elements of k are all identical to within a 
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Figure 3. The two nonzero elements of k^ for the left current 
are shown in panels A and B, with Tu — 0.5, Tfi = 1 and 
V — 4r. In each panel, the time dependence of the k^ element 
is shown at a range of band widths. 



Figure 4. The two nonzero elements of Ki, are shown in panels 
A and B, with Tu = 10, flc = lOF and V = AT. In each panel, 
the time dependence of the k^ element is shown at a range of 
inverse temperatures p. 



sign (koo.oo = Kii.ii = -Kii,oo = -koo,ii)- In Fig. [T] we 
therefore explore the dependence of one arbitrarily cho- 
sen element of k on a range of band parameters: cutoff 
energies Vic (the bandwidth is ~ 2Vlc) and band cutoff 
widths -. In each panel we go from a small bandwidth 
(red) to a large one (blue) at a set cutoff width, with the 
sharpest cutoff shown in panel A, an intermediate value 
in B and the smoothest in C. 

The effect of the the two parameters describing our 
chosen band shape on the memory kernel can be under- 
stood quite well by considering the trends shown in the 
plot: as V decreases, reflections are softened by the grad- 
ual slope at the band edge and the memory kernel decays 
more quickly. On the other hand, increasing the band- 
width induces oscillations at a frequency w « f^Ci but 
also increases the proportional weight of the short-time 
part of the memory kernel. Eventually, if we were to ap- 
proach the wide band limit, the memory would approach 
the form of a delta function and a Markovian description 
of the dynamics would become exact. 

The current memory kernel k^ is not as highly symmet- 
ric as K, and depends to some extent on all the parameters 
of the problem. There are two distinct (though similar) 
elements in ac^ at nonzero voltage, and these are plotted 
for the left current with Tv = 10 in panels A and B of 
Fig |2] The figure illustrates the dependence of k, on the 
bandwidth, which exhibits the same properties observed 
in K. This is also true of its ly dependence: to exemplify 
this. Fig |3] displays the same data for Ti/ = 0.5, where 
Kj decays more quickly and smoothly. Interestingly, the 
timescale over which k^ decays to zero does not appear 
to differ markedly from the corresponding timescale for 
K at similar parameters; this suggests that the cutoff ap- 
proximation remains as useful for the current as it is for 
impurity observables. 
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Figure 5. The two nonzero elements of k^ are shown in panels 
A and B, with Tv = 10, ^c = lOF and r/3 = 1. In each panel, 
the time dependence of the Ki element is shown at a range of 
bias voltages V . 



The effect of temperature on n^ depends greatly on 
the choice of other parameters. At the parameters we 
have chosen for Fig. [4] the asymmetry between the two 
Kt elements is increased somewhat when the temperature 
is lowered, corresponding to an increase in the current 
(not shown). One would expect a rather different effect 
when, for instance, things are set up in such a way that 
thermal enhancement of the current occurs. Unlike k, 
Ki depends on temperature even in the fully symmetric 
and noninteracting case is interesting, and expresses the 
fact that this quantity is connected to bath observables 
as well as to those in the impurity subspace. 

Fig. [5] shows how voltage affects the current memory 
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kernel: at zero voltage the two elements of k^ are identi- 
cal up to a sign, and the application of a voltage increases 
the diagonal element while suppressing the off-diagonal 
terms. Additionally, an increase in the oscillation fre- 
quency is observed for koo,ii, but no such clear trend 
exists for the oscillation frequency in kqo.oo- As V passes 
the bandwidth (~ 20r here), the left lead becomes en- 
tirely occupied and the right entirely empty, and further 
increasing the voltage ceases to have any effect on the 
current memory kernel, just as occurs in the case of the 
current itself (not shown). 

To show that the generalized NZME formalism in- 
troduced in this work indeed reproduces the correct 
results for the current as a function of time. Fig. [6] 
presents a comparison between the current obtained di- 
rectly (lighter solid lines) and by way of Eq. ( 30 1 (darker 
dashed lines). Pairs of lines describing the two differ- 
ent ways of obtaining currents at identical parameters 
overlap to within numerical errors, expressing the equiv- 
alence between the two methods when convergence in the 
cutoff time has been attained and the correctness of the 
approach at the tc ^ oo limit. 

Finally, while the NZME memory kernel technique and 
the cutoff approximation have been shown to be efh- 
cient for g i n a variety of interacting and noninteracting 
cases,HEll28l meaning that results at long times i 3> tc 
converge at a finite tc, no such calculations have pre- 
viously been carried out for the generalized technique. 
This entails a convergence analysis of the type exempli- 
fied graphically in Fig. |7] In essence, the cutoff time tc 
must be increased until the desired accuracy is reached. 
In the example shown in Fig. [7] convergence is achieved 
quickly and even short-time oscillations beyond the range 
of tc are predicted with some accuracy (as can be seen 
from the extension of the oscillatory ridges beyond the 
boundary of the transparent t — tc plane) . As an alter- 
native (if partial) representation of this idea, in Fig. ^ 
several plane cuts through this function are shown, but 
this time at Tv = 0.5; both the current and its time 
derivative are displayed. The rapid convergence visible 
in either representation illustrates that the idea of the re- 
duced dynamics technique and the cutoff approximation 
remains useful in practice even for the current, despite 
it being a non-system operator not accessible within the 
confines of the standard NZME formalism. 



VII. SUMMARY AND CONCLUSIONS 

We have reviewed the process of implementing reduced 
dynamics techniques by way of the NZME and the mem- 
ory cutoff approximation, which have recently been intro- 
duced with great success into several numerically exact 
nonequilibrium impurity solvers. The procedure of deriv- 
ing a memory kernel scheme for a general impurity model 
and obtaining calculable expressions in terms of standard 
second-quantization operators was outlined, and the ex- 
ample of the noninteracting resonant level model was 
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Q. =5. or mem 
Q =5. or direct 



rv=i 



— I 1 

a^=10.0rmem 
a =10.0r direct 



4- 



+ 



rv=l/2 ^^=15.0rmem 

---^=15.0r direct 
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Figure 6. The left current is shown as a function of the phys- 
ical time i at a variety of band parameters. Solid lines in 
Hght colors correspond to exact results calculated directly. 
Each such line is paired with a dashed line in a darker color 
at the same parameters showing converged results obtained 
from tracing over the t operator obtained from solving the 
generalized NZME Eq. (301. In all cases we have set F/? = 1 
and V = 4r. 




Figure 7. The time derivative of the left current in the cutoff 
approximation is shown as a function of both physical time t 
and the cutoff time tc. The transparent plane marks t — tc, 
and for t < tc (to the left of the plane) the results are exact. 
To converge the results for the current within some numerical 
accuracy, one must increase tc until —^ ceases to vary within 
that accuracy. Parameters are Tf = 10, ilc ~ lOF, T/3 = 1 
and V = 4r. 



worked out in full detail. For this noninteracting case, 
some illustrative examples of the physical properties of 
the memory kernel were discussed. 

An important limitation of the reduced dynamics tech- 
niques so far has been the lack of access to none-impurity 
observables, such as the electronic current in the resonant 
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Figure 8. The left current (top) and its time derivative 
(bottom) are shown as a function of the physical time t at 
Tv = 0.5, Qc = lOr, r/3 = 1 and V = 4r, for a range of cut- 
off times tc. The final line, labeled tc — oo, shows the exact 
result for comparison. 



bandwidth, voltage and temperature. The vahdity of the 
cutoff approximation for the current memory kernel was 
then verified and discussed. 

Looking forward, we expect the ideas expounded upon 
here to have several major implications: first, we hope to 
see them become a standard part of the toolbox of high 
quality time domain numerical simulations of impurity 
models, and ex tended to a variety of models and meth- 
ods. Second, since access to current enables access to 
Green's functions, the benefits offered by memory tech- 
niques are expected to be applicable to interacting lattice 
simulations as well, by way of mapping schemes such as 
dynamical mean field theory and its various extensions. 
Finally, we believe the memory kernel framework is a fer- 
tile ground for defining new approximation schemes more 
general than the cutoff approximation, and in this con- 
text it will be particularly interesting to understand the 
long-time behavior of the memory kernel in interacting 
cases and its behavior in larger impurity models. 
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